Spherically symmetric spacetimes in f{R) gravity theories 
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We study both analytically and numerically the gravitational fields of stars in f{R) gravity theo- 
ries. We derive the generalized Tolman-Oppenheimer-Volkov equations for these theories and show 
that in metric f{R) models the Parameterized Post-Newtonian parameter 7ppn = 1/2 is a robust 
outcome for a large class of boundary conditions set at the center of the star. This result is also 
unchanged by introduction of dark matter in the Solar System. We find also a class of solutions 
with 7PPN « 1 in the metric f{R) = R — fJ-"^ /R model, but these solutions turn out to be unstable 
and decay in time. On the other hand, the Palatini version of the theory is found to satisfy the 
Solar System constraints. We also consider compact stars in the Palatini formalism, and show that 
these models are not inconsistent with polytropic equations of state. Finally, we comment on the 
equivalence between f{R) gravity and scalar-tensor theories and show that many interesting Palatini 
f{R) gravity models can not be understood as a limiting case of a Jordan-Brans-Dicke theory with 
uj -» -3/2. 

PACS numbers: 04.50.-|-h, 98.80.-k, 95.35.-|-d 



I. INTRODUCTION 

The observation that the expansion rate of the universe 
appears to be accelerating [1, 2] has led to a great inter- 
est in exploring the possible extensions of the Einstein- 
Hilbert theory of gravity. A particularly popular class 
of models involves including nonlinear interactions in the 
Ricci scalar R: 
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where 6*111 is the matter action. Setting f{R) = i? — 2A 
corresponds to the canonical Einstein-Hilbert action with 
a cosmological constant A. Recent interest in these mod- 
els has followed from the observation 0, |3| that adding a 
function Sf{R) = f{R) — R = —jj^/R can give rise to the 
observed acceleration without a cosmological constant. 

However, it has been proven that when understood as 
a metric theory [3l|, the action H]) leads to predictions 
that are in contrast with the measurements in the Solar 
System d, Q. Indeed, as was claimed by Chiba [sj, the 
Parameterized Post-Newtonian (PPN) parameter 7ppN 
in these theories is 7ppn — l/^: while the observational 
constraint requires 7ppn — 1 ^ 10~^J3L This result 
has been contested by several authors [aH, , but was 
recently confirmed by a more direct computation by Er- 
ickcek et al. Q- Both Chiba and Erickcek et al. consid- 
ered the space-times of stars surrounded by empty vac- 
uum where at large radii R V^fJ,'^ = i?vac for the 
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-fjr/R model. However, in reality stars and 



'Electronic address: Kimmo.Kainulaincn@ phys.jyu.fi] 
^Electronic address: Johanna. Piilonen@phys.jy u.fi| 
Electronic address: Vappu. Reijonen@helsinki.fiT 
§ Electronic address: |Daniel.Sunhede@phys.jyu.fi| 



star-systems in galaxies are surrounded by a halo of dark 
matter and one could ask if this might change the con- 
clusions. 

Indeed, in General Relativity (GR) one finds R ~ 
SttGp (assuming p <^ p) and this relation becomes the 
appropriate asymptotic limit for the Ricci scalar in many 
,f{R) models with a constant density background. Typ- 
ically one assumes that pdm ~ 0.3 GeV/cm'^ [ll|, so 
that SirGpuM ~ 10®i?vac- This continuous density field 
would then appear to have a potentially large influence 
on the predictions of an f{R) model [S^l: for the case 
f{R) = i? - p-^/R one has 
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and according to the previous argument one expects that 
p'^/R^ ~ P^/8ttGp ^1. It would then be tempting to 
conclude that the theory reduces to GR in a high den- 
sity environment. If so, then 7ppN ~ 1 in accordance 
with GR but in contrast to the vacuum prediction of the 
theory. This argument was essentially what was used in 
a recent attempt [l3] to show that also metric theories 
predict 7ppN = 1 in a dark matter background, in con- 
trast with the vacuum result 0, Q . We show here that 
introducing a finite dark matter density around a star in 
fact has no effect on the solutions in metric f{R) gravity. 

However, another issue with the metric f{R) theories, 
which contain derivatives of fourth order, concerns the 
boundary conditions. Typically the boundary conditions 
are set far outside the star for the exterior solutions and 
to our knowledge no interior solutions have previously 
been computed in the literature. Our approach here is 
to set the boundary conditions at the center of the star 
and compute the solutions throughout, both interior and 
exterior to the star. We find that the result 7ppn ~ 1/2 is 
robust for a large class of boundary conditions, at least as 
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long as one requires that the metric is finite at the center 
of the star. These solutions are also unaffected by intro- 
duction of dark matter in the Solar System. We also find, 
however, that another class of more fine tuned solutions 
exist for which the Ricci scalar settles in the Palatini 
limit of the trace equation (defined as when derivative 
terms become negligible), giving 7ppN ~ 1 even in the 
metric theories. Unfortunately these solutions are not 
physically acceptable since they are unstable in time, de- 
caying through the tachyonic instability discovered by 
Dolgov and Kawasaki [ij]. We will argue that our re- 
sults are not restricted to the f{R) = R ~ fJ.'^/R model, 
and hence generic metric /(i?) theories appear to be ruled 
out by the PPN data. 

We will also consider finite density effects on stellar 
solutions in the Palatini formalism. In this case the limit 
R « SnGpBM does hold, reducing the theory to GR not 
only near and inside stellar objects, but essentially at 
all scales where the matter density is significantly higher 
than the asymptotic value of the cosmological constant. 
As a result. Palatini f{R) gravity leads to 7ppN = 1 to 
a very high accuracy in the Solar System. These solu- 
tions are also stable in time since the Dolgov-Kawasaki 
instability does not exist in the Palatini formalism p^ . 

Barausse et al. [15] recently claimed that no solu- 
tion can be found for Palatini f{R) models consistent 
with a polytropic equation of state (EOS) with index 
3/2 < F < 2. We do not agree with the interpretation 
of their results. The problem discovered in Ref. [l^ fol- 
lows from the fact that the Einstein equations in f{R) 
models depend on p" through the r-derivatives of the 
function F = df/dR, and that this quantity diverges at 
the boundary of the star for 3/2 < F < 2. However, 
we will show that this singularity is so weak that, in the 
case of a neutron star, it becomes relevant only at dis- 
tances Ar < 0.3 fermi from the boundary. Defining a 
sharp boundary to this degree of precision is obviously 
unphysical and beyond the validity of any polytropic 
model. Moreover, a polytropic description will always 
brake down at scales less than the mean collisional dis- 
tance of the particles in the fiuid so there is in fact no 
physical singularity. As a result, at least all f{R) models 
where F = df /dR is a decreasing function of R should be 
consistent with the existence of compact GR-like stars. 

Finally, we will consider the scalar-tensor gravity the- 
ory (STG) equivalence for metric and Palatini f{R) grav- 
ity. The Palatini formalism in particular corresponds to 
a STG model with no kinetic term in the Einstein frame, 
and as a result the Palatini field follows the shifting po- 
tential extremum with no resistance. This explains why 
Palatini models can be GR-like, while their metric coun- 
terparts fail to obey the PPN constraints. However, Pala- 
tini f{R) models whose effective potential in the Einstein 
frame has negative curvature, V^g < 0, are found to be 
problematic in the sense that they can not be under- 
stood as a limiting case of Jordan-Brans-Dicke theories 
with uj —* —3/2. 

The paper is organized as follows: sectionlTTlintroduces 



a complete set of Tolman-Oppenheimer-Volkov equations 
for solving the configuration of a static spherically sym- 
metric star in metric f{R) gravity. We derive the New- 
tonian limit of these equations and compute analyti- 
cally the predicted value for 7ppn- We also supplement 
the analytical work with a direct numerical solution us- 
ing the complete Tolman-Oppenheimer-Volkov equations 
with varying boundary conditions and different forms of 
f{R). In section IIIII we consider the Palatini version 
of the theory and show by both analytic and numer- 
ical calculation that 7ppn — 1 for any Palatini f{R) 
model where F = df/dR is a decreasing function of R. 
We also consider compact Palatini stars with polytropic 
equations of state in this section. In section [TVl we con- 
sider the scalar-tensor equivalence of f{R) theories and 
section IVl contains our conclusions. 

II. METRIC f{R) GRAVITY 

In the metric formalism one assumes that the Ricci 
scalar R and the covariant derivative are given in 
terms of the Levi-Civita connection. The field equations 
are thus obtained by varying the action ^ with respect 
to the metric g^j/ only: 

^^M- - lf9^^^ - V^V.F -f g^.UF = 8^GT^, , (3) 

where F = df/dR and □ = g^^V fjVi,- In the spheri- 
cally symmetric case of a static mass distribution, these 
equations will contain only three independent parameters 
corresponding to the Ricci curvature R and the tt- and 
rr-components of the spherically symmetric metric g^^: 

ds^ EE g^^xf^x" ^ -e^'^'-^dt^ + e^('')dr2 + r^dfl^ . (4) 

To solve the functions R, A and i?, it is most convenient 
to use the tt- and rr-components of the field equations 
supplemented by the trace equation for R: 

DF+^{FR~2f)^^T, (5) 

where T = = —p + 3p. In the static and spherically 
symmetric case Eqn. ([5]) becomes 

2 A' - R' 

F" + -F' + - —F' + —{FR-2f)^—%nGT. (6) 

T" 2 3 3 

Note that if the derivatives of F (eventually of R) vanish, 
Eqn. (O reduces to 

Fi? - 2/ = SttGT , (7) 

which for f{R) = R — /i^/i? has the solution 

Rt = ^ (-SvrGT ± V(87rGT)2 + Up^^ . (8) 

Given that p <^ p and /x^ <C SnGp, the solution with 
a positive sign, corresponding to an asymptotically de 
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Sitter space gives R « SnGp as mentioned in the 
introduction. This argument suggests that, in analog to 
the vacuum hmit of Ref. [1], the solutions for Eqn. ([6]) 
should smoothly approach the solutions of ([7]) for large 
enough r in a finite, spatially constant background den- 
sity. 

It is useful to rewrite the field equations ([3]) in the 
equivalent form: 
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where G^i, = — ^Rg^^v The source equations for A' 
and B' can now be found from the tt and rr components 
of the field equations ©. Using Eqn. ^ to get rid of 
the F" term on the right hand side in the rr equation, 
one eventually finds: 
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+ ^(^+|;)-7^', (11) 

where 7 = rF' /2F . Furthermore, from the conservation 
equation V^T^'' = (valid in the Jordan frame one 
finds the following simple relation between p' and A': 
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P^-ylP + P)- (12) 
When supplemented by an equation of state 
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equations ([6]), (|10mip and give a complete general- 
ization of the Tolman-Oppenheimer-Volkov equations for 
a relativistic, spherically symmetric star in metric f{R) 
gravity theory. The above equations, which promote F 
and F' to be free variables along with A and B, appear to 
be much simpler than the alternative equations derived 
e.g. in Ref. [18], where R and eventually F and F' were 
expressed in terms of A and B and their derivatives. 



A. The Newtonian approximation 

As was emphasized by Erickcek et al. a unique 
exterior solution for a stellar object is found by match- 
ing it with an interior solution in the presence of matter 
sources. In this section we find such a solution for the 
Tolman-Oppcnhcimer-Volkov equations derived in the 
previous section. If any solution is to be found compati- 
ble with the PPN limits, the corresponding metric should 



be close to that of the Newtonian limit. That is, one can 
assume that ^, B < 1 and A', _B' < l/r in Eqns. ^ and 

(Hmni). 

Another related simplification that we shall make is 
to neglect the pressure throughout and use fixed den- 
sity profiles instead. The rationale for this is that we 
are ultimately only interested in the metric created by a 
generic Newtonian matter configuration and not on solv- 
ing for these configurations self-consistently. Turning the 
argument around, Eqns. (|12m3p could be used in retro- 
spect to compute the pressure and equation of state that 
would be needed to create the employed density profile, 
but since p <^ p in the Newtonian limit, that pressure 
would only give a negligible contribution to the metric 
and can hence be ignored. In fact, except for imple- 
menting /(i?) gravity into a full numerical stellar model, 
our approach may provide the most accurate test for the 
model, since the equation of state computed in retrospect 
from a realistic stellar density profile is likely to be more 
appropriate than simple toy model equations of state (see 
section (|III A|) below). 



1. The trace equation 

Let us first consider the trace equation ([6]) in the New- 
tonian limit: 

F" + lF'+^-{FR~2f)^-^p, (14) 

where we have also assumed that p <^ p so that T ^ —p. 
Furthermore, define a new variable d = F — 1 which 
parametrizes the deviation from GR. For the particular 
model under consideration, f{R) = R — p'^/R, one now 
obtains the following equation for d: 
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Note that, save for the nonlinear term, this form is of 
course independent of the chosen f{R) model. That is, 
the class of the solutions for d (and hence for F and 
ultimately for R) is completely dictated by the relative 
size of the nonlinear term to the derivative terms. It 
is hence instructive to rewrite the nonlinear term and 
the source in (jlSp so that their relative size in different 
environments become more transparent: 



2^, fm^{i-3d) 



\ SVd 



P 

Pcrit 



(16) 



where 5Hq = SnGpcrit is the Hubble expansion and 
V3/i^ = 4A. Assume now that one starts to integrate 
Eqn. (fTB]) outwards from the center of a stellar object 
with some boundary value d = do. Obviously, the den- 
sity ratio pq /pcrit is enormous at r = and it completely 
dominates the evolution unless do is fine tuned to zero to 
a high precision. Indeed, the nonlinear term can influ- 
ence the evolution of d only at very low densities, where 
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P Pa = f^APciit, or at high densities when d < dp, 
where 



(17) 



and Rp is given by Eqn. ([8]) under the assumption that 
pressure is neghgible. At the center of the Sun one has 
Po/Pcrit ^ 10'^^, and thus the nonhnearity is neghgible 
unless d < 10~^^! At the edge of the Sun (defined as 
the radius that encompasses 99.9% of the mass) the cor- 
responding value would be d < 10~^^. Even outside 
the Sun, within the local dark matter distribution with 
Pdm = 0.3 GeV/cm'^, one finds that the density will dom- 
inate over the nonlinear term unless d < 10~^^. 

Barring such fine tuned boundary conditions, the non- 
linear term can be neglected in Eqn. psp . at least during 
the initial stages of the evolution. Let us now solve for 
d in this approximation, setting the boundary conditions 
at r = ro with the understanding that in the end tq —> 0. 
Integrating Eqn. once, one finds: 
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where ( )o indicates that an expression is evaluated at 
r = To, Am(r) = m(r) — m(ro), and 



m(r) — 



P- 



(19) 



Eqn. pS)) shows that the boundary value for d' at 
To = plays no role in the solution as long as it is finite, 
because then (?'^c?)q at ro 0. Only considering 
finite solutions, one can set ro = and (r^d')o = and 
perform another integration to get: 
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Another important observation to be made from 
Eqns. ([T5)) and PPI) is that d gains a negative contribution 
from the density evolution. Thus, even if one starts from 
a value d^, for which the nonlinear term is negligible, 
the evolution will push d towards the region where non- 
linearity does become important. Whether this happens 
inside the star or a stellar system depends on the relative 
size of the matter effect compared to do, and can easily 
be estimated in retrospect. For now, assume that do is 
large enough so that nonlinear terms can be neglected 
throughout. 

Outside a star in the Solar System scale, the total mass 
of the configuration is overwhelmingly dominated by the 
mass of the star: given pdm = 0.3 GeV/cm"^, the mass of 
the dark matter component inside a radius r in the Solar 
System is roughly Mdm ~ IO'^'^Mq (r/AU)^, where AU 
is one astronomical unit. This shows that the exterior 
finite dark matter density has no effect on the solutions, 
contrary to the claims made in Ref. ^12i] . Hence, to an 
excellent approximation, the exterior solution is 
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where M is the total mass of the star, M = m{rQ). The 
constant has been C extracted for convenience so that 
the 1/r law becomes explicit in (|2T|) . 
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Putting in the numbers for a typical star like our Sun, one 
finds that the scale of the density induced evolution in d is 
of order d—do ~ 10~^ — 10~^, from the center to distances 
on the order of AU. This sets the scale for the validity of 
the approximation of neglecting the nonlinear terms: we 
conclude that the solution i21\) holds for any dg > 10~^ 
and for arbitrary finite d'^ at ro = 0. Furthermore, it is 
not affected by a finite dark matter distribution external 
to the star. However, for boundary conditions c?o ^ 10^^ 
a more careful analysis is needed and we will study these 
solutions in separate in section fll CI 

We have based our argumentation on the idea of inte- 
grating the equations starting from the center of the star 
with some boundary condition. The Erickcek et al. solu- 
tion 6], on the other hand, was derived with a boundary 
condition R —> \/3p^ or equivalently d — > 1/3 far away 
from the star at r — > cx). However, it is obvious from the 
previous discussion that this boundary condition corre- 
sponds to an almost identical boundary condition at the 
center of the star: do = 1/3 — C w 1/3, up to correction 
of order 10~^. This also appears to be the physically 
best motivated boundary condition, as stars were born by 
condensation of inhomogeneities from an initially rather 
smooth matter distribution, whereby one most naturally 
expects that R ~ a/S/x^ at large scales. 



2. Source equations for A and B 

Let us now turn to the source equations for A and B in 
the Newtonian limit, armed with the solutions ((20)) and 
dSI]) for F = 1 + d. For this solution, 7 = rF' /2F « 
rd'/2 = —GM/3r ^ 1, so that one can expand in this 
variable as well. Hence, to first order in small quantities: 



A' 



{rBY 



B 



2d' 



II r 



r 

IQttG 



r'^p + 



-d ' 
(Mr)2 



3Vd{l + d) 



(23) 
(24) 



where we have neglected pressure but kept the nonlinear 
contribution, rewritten in terms of d = /Lt^/i?^. However, 
given the solution for d, it is easy to see that the 
nonlinear terms are always completely negligible in com- 
parison with the other terms, at least as long as one starts 
with do ^ 10~^. Neglecting these terms, it is straightfor- 
ward to obtain the solution 
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Bo+Ao, (26) 
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where we have kept all integration constants, setting the 
boundary values at r = tq. In particular, one finds that 
the exterior solution is: 
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where C is given by Eqn. dm [11] ■ Note that the validity 
of the Newtonian limit for the above solutions can be 
verified in retrospect from the final forms of the solutions 
(I^Tj) and (|27ll28p . Apart from arbitrary constants: A, B, 
d^GM/r<$: 1. 

In order to obtain the appropriate asymptotic limit 
A ^ 1/r for large r, one must set Aq = —AC — i?o in 
Given this constraint one finds in the Newtonian 

limit: 
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Thus, for any finite value of Bq at the center of the star 
one obtains 



1 



7PPN 



(30) 



This is of course the well known result by Chiba , which 
was later confirmed in Ref. 3 and yet again contested 
in Ref. We have shown here that this result holds 

for a large class of boundary conditions at the center of 
the star: when dg > 10~^ and for practically any dg, Aq 
and Bq, as long as one remains in the Newtonian limit 
and A ~ 1/r for large r. Furthermore, this solution is 
not affected by the dark matter distribution surrounding 
the star. 

Turning the above argument around, one can also sec 
from Eqn. that the condition 7ppN = 1 exterior to 
the star necessarily leads to a blow-up in the metric com- 
ponents at the center (again assuming that do ^ 10^'^); 
indeed, 7ppN ~ 1 would imply —Aq ~ -Bq ^ GM/rg 
CO when ro 0. Such a solution is clearly unphysical, 
even though the actual divergence of the metric occurs 
beyond the validity of the Newtonian approximation. 



B. Numerical analysis 




FIG. 1: Shown is the density profile of the star (solid) and 
the appropriate pressure needed to support this configuration 
(dashed) (in units ?"q^). The inset figure shows the equation 
of state p = p(p) corresponding to this solution. 



figure shows the pressure profile computed in retrospect 
from equations (|12p and (|13p . using a known solution for 
the metric. The inset figure displays the corresponding 
equation of state p = p{p). These plots show that the as- 
sumption p <^ p holds, and that our profile corresponds 
to a reasonable equation of state. 

Fig. ^ shows the evolution of A, B, and d — dvac, 
where dvac is the asymptotic value of d in vacuum. These 
plots turn out to be practically independent of the choice 
of boundary conditions do, d^, Bq and Aq, as suggested 
by our Newtonian analysis. (One of course has to set 
Aq — —Bq — AC in order to ensure the appropriate limit 
for A at large r). The evolution of the corresponding 
7PPN parameter is shown in Fig. ([3]). Although it of 
course only makes sense to talk about 7ppn for the exte- 
rior part of the solution, it is nevertheless interesting to 
compare how the quantity evolves in GR and in a met- 
ric f{R) model, reaching different constant values at the 
edge of the star. We find that this form of the solution is 
generic and independent of the boundary conditions as 
long as do ^ 10~^ and the metric is finite at the center 
of the star. 



We will in this subsection support our analytic re- 
sults (which were based on the Newtonian approxima- 
tion) with numerical solutions of the exact source equa- 
tions (jlOmip and the trace equation ([6]), still using a 
fixed density profile. The actual profile used corresponds 
to the known density profile of the Sun with a central 
density of 150 g/cm^ and with a roughly exponential de- 
pendence on r. The dark matter distribution was taken 
to be a constant with pdm = 0.3 GeV/cm'^, superim- 
posed on the profile of the star. We display this density 
profile in Fig. ([1]) (solid line). The dashed line in the 



1. Other forms of f(R) 

We have repeated the above numerical analysis for sev- 
eral other models in addition to f{R) = R — p'^/R. As 
suggested by the Newtonian analysis, all f{R) models 
should give 7ppn — 1/2, as long as the nonlinear terms 
can be neglected. We have verified this result by nu- 
merical analysis, where have also considered explicitly 
various other models, including f{R) = R — fj,'^/R + aR^ 
and f{R) = R — PR"' , and in particular a model similar 
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Q 2 4 6 8 10 corresponding solution in GR (dashed). 

r/r0 



FIG. 2: Shown are the functions A (red) and B (blue) for the 
metric f{R) = R — /i^/-R model (solid) and for GR (dotted). 
Also shown is the function d — dvac (dashed green) , where dvac 
is the asymptotic value of d in vacuum. 

to the one suggested in Ref. [l^: 

f{R)^R + aVR, (31) 

as a possible candidate for passing the Solar System con- 
straints. We find that all models fail the PPN limit; as 
long as the model parameters are set to give the correct 
asymptotic cosmological constant and one does not add a 
true cosmological constant to the f{R) function, all mod- 
els produce results that are essentially indistinguishable 
from the f(R) = R — jJ."^ /R model in the Solar System 
scale [ssj . 

C. Solutions with do < IQ'^ 

Let us now go back to the class of solutions with very 
small boundary values for do. As explained above, in this 
regime the matter induced evolution of d is strong enough 
to push the solution to the nonlinear region inside a Sun- 
like star. One can argue qualitatively that the resulting 
solution will be one where d oscillates around the value of 
dp corresponding to the Palatini limit, Eqn. ([5]). Indeed, 
since i? is a real number, d must always remain positive 
in the f{R) = R — n'^/R model. However, starting from 
rfo S> dp at r = rg, d will first start to decrease. This 
evolution is bound to be reversed by the nonlinear term 
before d becomes negative, but once d starts to increase, 
the nonlinearity shuts off again and the the finite density 
effect turns the evolution back towards smaller d. As the 
cycle gets repeated, the result is an oscillatory motion 
around the Palatini limit, defined as the solution of the 
equation ([7]). Wc show an example of this behaviour in 
Fig. Q. The solid line shows the evolution of d, which 
indeed settles to a damping oscillatory pattern around 
the Palatini limit, shown by the dotted line. We also 



display the metric coefficient A (dashed) which, after a 
short interval of "Newtonian f{R) evolution" , settles to a 
converging oscillatory track around a path parallel to the 
GR solution ^gr — ~ 0. The solution for B turns out 
to be numerically indistinguishable from the correspond- 
ing GR solution. In summary, A and B turn out to be 
very close to the GR solution simply due to the fact that 
the Palatini solution is virtually indistinguishable from 
the GR metric (see section llll| . 

Note that at the center of the Sun, the oscillations oc- 
cur in scale ~ lO~^*r0, so it is not numerically feasible to 
continue the solution all the way to the surface. We have 
nevertheless run the code over thousands of oscillation 
periods, verifing that the solution does indeed stabilize 
around the Palatini limit. Furthermore, this behaviour 
is independent of the boundary value do. (Of course, if 
one sets do = dp, the solution will become flat without 
any oscillations.) The above example used a very small 
do, but the qualitative behaviour of the solution should 
remain the same for any do for which the Newtonian evo- 
lution is strong enough to bring d to zero inside the star. 
As a result, it is safe to conclude that for sufficiently 
small do the solution will be such that inside and in par- 
ticular outside the star A « >1gr and B « -Bgr, so that 
7ppivj « 1. In practice, the boundary for this result may 
be somewhat less than d < 10~^ since d needs to reach 
the nonlinear region already close to the center of the 
star. If not, the initial evolution of A and B will have 
time to push the metric and eventually 7ppn too far from 
the GR solution. 



1. The Dolgov-Kawasaki instability 

The above section explored an attractor solution 
around the Palatini limit for small values of the bound- 
ary value do. However, it turns out that this class of 
solutions is related to the well known Dolgov-Kawasaki 
instability in the f{R) = R — fJ,^/R model. Perturb- 
ing around the static solution, d(r) d{r) + Sd{t,r), and 
expanding to first order in the perturbation one obtains 
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FIG. 4: Shown are the functions d (sohd) and A—Aq (dashed), 
where Ao is chosen such that for GR A^l/ratr— »cx3. The 
dotted line shows the Palatini value of d, which corresponds 
to STvGp. 



the following equation: 



d^Sd - V^Sd 



6d3/2 



{l + 3d)Sd = -mlSd, (32) 



The negative effective mass squared in this wave equation 
is what, when restricted to the static limit, gives rise to 
the damping oscillations in r around the Palatini limit 
seen in Fig. ([4]). However, this attractor behaviour in 
r comes with the price of making the solution unstable 
in time. Indeed, expanding the perturbation in Fourier 
modes, one finds that a mode with wave vector k has 
time dependence (for < 0) 



Sdk{k, t) 



(33) 



so that there are unstable modes with k < |mrf|. This is 
the instability first found by Dolgov and Kawasaki ^13,] . 

It is important to notice that the magnitude of d con- 
trols both the characteristic time scale of the instability, 
iinst ~ and the shortest scale of the unstable 

modes, Hnst c^inst- For the GR-like solutions above 
4 X 10~^^ at the center of the Sun and one 



8 X 10" 



sec. Even for 



finds an unfavourable tji 
^ ~ dpDM ~ 3 X 10~^^ corresponding to the dark mat- 
ter density, the time scale of the decay is still relatively 
short, tinst ~ 6 yr. The Dolgov-Kawasaki instability then 
clearly forbids the Palatini tracking solutions as physical 
ones. For solutions with d ^ do > 10^^ the instability 
time is ii„st > 500 kyr, and it becomes of order 1 Gyr 
for d ~ I. In summary, one can conclude that the metric 
f{R) = i? — /i^/i? model has both approximately stable 
solutions with 7ppN = 1/2 and unstable solutions with 
7PPN — 1 , but that it does not have any sufficiently stable 
solutions that would also pass the Solar System tests. 

Quickly after the discovery of the Dolgov-Kawasaki in- 
stability, a way to cure it was suggested in Ref . ^21i] . The 



idea is to add for example a quadratic term (a/2/i^)i?^ to 
the theory, after which the effective mass in equation 
([5^ becomes (for d < 1): 



3a - 6d3/2 



(34) 



Since the theory is stable if is positive, could one 
stabilize the above GR-like solutions in this way? A pos- 
itive thing to this end is that the GR-like solution does 
remain an attractor for boundary values do < 10~^ even 
in this extended model. However, in order to make the 
theory acceptable, one has to make sure that the mass is 
positive everywhere in the Solar System. That is, 
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a > 2^3/2 



IQ- 



(35) 



Unfortunately, even for this small value of a, the function 
F, which controls the density influence on the growth of 
A and B in Eqns. (|10j|lip . becomes enormous inside the 
star: 



1 



aRp 
1^ 



SnGp 

^ 2~ 



lo"-lol^ 



(36) 



for p ^ pq 0.1 — 150 g/cm'^. It is obvious that such 
a value of F would completely shut off the evolution of 
A and B, giving rise to a nearly massless Schwarzschild 
exterior solution. This argument is quite generic and 
it would thus seem to be difficult, if not impossible, to 
stabilize a GR-like solution throughout a stellar system 
in any metric f{R) model. Note that a ~ 1 stabilizes the 
model at all scales. 



III. PALATINI f{R) GRAVITY 

We have recently studied Solar System constraints on 
Palatini f{R) gravity in Ref. |lB|, assuming that the star 
was surrounded by vacuum. This section will include the 
effects of the nonzero dark matter density to the analysis. 
As is well known, the trace equation in the Palatini case 
is an exact algebraic equation identical to ([7]): 

FR-2f = SttGT , (37) 

where the Ricci scalar is a function of both the metric 
and the independent affine connection, R — g^" Rf^,y{r). 

Td^ . the source equations for 



As was was first shown in 
A{r) and B{r) become: 



A' 



-1 



B' = 



1 + 7 
1 

1+7 



—SnGrp 

r 

+ —SnGrp 
F 



where 




(38) 
(39) 

(40) 
(41) 
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and 7 = rF' /2F as before. The conservation equation 

relating p' and A' and the equation of state 
remain unchanged. 

Once again, if any solution is to be found compatible 
with the PPN limits, the corresponding metric should 
be close to that of the Newtonian limit. That is, one 
can assume that A,B <^ 1 and 7 <C 1 in Eqns. (|38t - 
159]) . From the trace equation (|57|) one knows that for 
the f{R) — R ~ ^J''^ /R model, R ~ 8ttG so that also 
in the Palatini formalism d = F ~ 1 <ti 1 inside any 
matter distribution for which SttGT ^ 12 fp. Neglecting 
pressure, it is then straightforward to show that to first 
order in small quantities: 



A' 
(rB)' 



B 



2d' . 



SnGr^p+ir^d'Y 



(42) 
(43) 



Except for the terms containing d, these equations are 
identical to the corresponding equations in GR. That is, 
if the terms containing d are negligible, these models will 
indeed be compatible with Solar System constraints. 
The solution to the above equations is easy to find: 



B 
A 



2Gm{r) 



' +^^+rd'-^^. (44) 
TO r r 

2G r\r'm~^-^^B.^A. 



(45) 



where we have kept all integration constants setting the 
boundary values at r = tq. Now, for the specific model 
d — / R^ it is straightforward to show that for any 
density profile that is non-singular at the origin, (r^(i')o 



will go to zero at r : 
solution is: 



0. Hence, one finds that the exterior 



B 
A 



2GM (rB] 







rd' . 



r 

2GM 



r 

(.rB)„ 



(46) 



Bo + 3C + Ao + d, (47) 



where C is given by Eqn. (|22p . Furthermore, since 
d ~ ji^ /{S-KGpf and d' ~ [p' / p)d, these terms are in- 
deed completely neghgible compared to GM/r. That is, 
setting Ao = — 3C— i?o in order to obtain the appropriate 
asymptotic A ^ \/r for large r, one finally obtains: 



7PPN 



B 
A 



2GM + 3(rB)o 



2GM + 3{rB)o ' 
so that the model is indistinguishable from GR: 



7PPN 



1 



(48) 



(49) 



The above results are in no way surprising. As was 
shown in [l^ (see also [22|), for the particular model 
in question F is extremely close to 1 inside any density 




FIG. 5: Shown are the metric functions A and B for the 
Palatini f{R) — R ~ P*' /R model (solid), which completely 
overlap with the GR solution (dashed). 



distribution much larger than the effective cosmological 
constant in vacuum. Hence, d = — 1 is completely neg- 
ligible and the very same supression mechanism brings F 
to 1 also in the exterior, dark matter dominated region. 
That is, in the Palatini formalism, a star surrounded by 
dark matter will be even more GR-like than a star sur- 
rounded by vacuum. Observationally, both scenarios pre- 
dict a post-Newtonian parameter 7ppN indistinguishable 
from 1. 

To support the analytic results above we have also 
solved the exact source equations p8ll39p using the same 
approach and density profile as in the metric section. 
Fig. [5] shows the numerical solution for A and B in the 
Palatini formalism, which completely overlap with the 
GR solution. Hence, the curve for 7ppN in the Palatini 
formalism is identical to the GR curve shown in Fig. [3l 
As suggested by the Newtonian analysis, these plots turn 
out to be generic and independent of the given boundary 
conditions Bq and Aq. 

We have also repeated the numerical analysis for sev- 
eral other forms of f{R). We find that for all models 
where F = df/dR is a decreasing function of R, 7ppN 
remains indistinguishable from 1 as long as the model 
parameters are chosen to reproduce the correct asymp- 
totical cosmological constant. 



A. Polytropic stars in Palatini f{R) gravity 

A recent paper by Barausse et al. considered static, 
spherical stars with polytropic equations of state, and 
concluded that Palatini f{R) models would not give 
physically acceptable alternatives to GR. While the anal- 
ysis does reveal a genuine qualitative difference between 
Palatini f{R) gravity and GR, we do no t ag ree with their 
conclusions. The main criticism of Ref. llSlI was based on 
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the argument that for a star with a polytropic equation 
of state, 



r , P 

P = , p = po + 



r - 1 



(50) 



the polytropes with 3/2 < F < 2 lead to blow-up of F" 
and hence also of the metric at the "edge" of the star, 
f = '"out, where po — s- 0. It is indeed true that due to 
the F' and F" terms on the r.h.s. of Eqns. p8ll39p . the 
metric may be more sensitive to rapid density changes 
in f{R) theories than it is in GR. Since F" shows up 
in the equation for B' , F" blowing up at the boundary 
leads to a divergence of B' as well. However, we wish to 
show that this issue has more to do with the peculiarity 
of polytropic EOS's, when assumed to hold to a mathe- 
matically abstract precision at the boundary, than with 
the theory of Palatini f{R) gravity. 
Let us begin by writing 



F' 
F" 



OF , 
dp 
d^F 



9p2 



l\2 



OF 



(51) 



The partial derivatives OF /dp and d^F/dp^ can be com- 
puted by using the trace equation (|37p and are in gen- 
eral just some finite numbers at p = 0. The behaviour 
of F' and F" is then completely controlled by the be- 
haviour of p at the boundary. Since the continuity equa- 
tion is valid for any f{R) gravity theory in the Jor- 
dan frame, one finds that (given a finite and nonzero A') 
p' = {dp/dp)p' ^ {dp/dp){p + p), so that for a polytropic 
EOS (EOl): 



P 
P" 



2-r 



3-2r 



(52) 



in the limit p Q, independently of the gravity theory 
in consideration. In particular, for the polytropes with 
3/2 < r < 2, the second derivative p" diverges while p 
and p' remain continuous in the vacuum limit. Note that 
the behaviour in Eqn. ([5^ (for F > 1) is only possibleif 
the star has a sharp edge where the density goes to zero 
at some finite radius rout- 

It is easy to solve Eqn. ([5^ for the density profile close 
to the edge: 



p{r) = 



(rout - r)r-i , r< Tout 
, r> Tout 

N(3-2r)/(r-i)/ 



(53) 



r)y- ^'ViTont - r), 

Tout for F > 3/2. At 



One hence finds p" ^ (rout ^ 
which indeed blows up at r 
first sight, this looks troubling since certain interest- 
ing microphysical models predict polytropic EOS's with 
F G [3/2,2]. For example, a degenerate and nonrela- 
tivistic fermion gas has F = 5/3. However, for this to 
be a problem, the hard polytropic EOS would have to 
be valid all the way until p — 0, which is an unphysical 
expectation. 



First, more realistic EOS's tend to become softer as 
the density decreases. For example figure 5 in Ref. [23| 
shows that the effective adiabatic index F, which can be 
as large as 2.7 at the interior of a neutron star, drops 
to about 1.25 when p < 10^*^ g/cm'^. In such a case the 
metric will be perfectly well-behaved everywhere (36| . 

Second, since the function F" is typically controlled 
by a small parameter (for the f{R) = R — p^ /R model, 
F" ~ ^*), one would expect that the singularity at the 
edge should be very weak. To see this quantitatively, let 
us consider the tidal acceleration close to r = rout- For 
a static, spherically symmetric metric one can show that 
the equation of geodesic deviation in the radial direction 
is given by: 
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-{2A" + A'^ - A'B')ri'' 



(54) 



where r]^ is the radial spatial frame component of the or- 
thogonal connecting vector, D^/Dt^ = (w^V^)(f;''V^), 
T is the proper time, and v'^ is the tangent vector at each 
point of the timelike geodesic [33|. The angular compo- 
nents are proportional to e~^A' and will hence remain 
unaffected by the blow-up in F" . Using Eqns. (|38ll4ip 
and only keeping the terms proportional to F" , one finds 
the effect of the singularity on the tidal acceleration: 



K ■ 



-B j,2pll 
i V 



3e^A 



5 F' 3 /F' 



2rF 



F 



(55) 



where X = ^{R — f / F) is the radially dependent effective 
cosmological constant. 

A reasonable approximation for the magnitude of A' 
close to the surface is A' ^ 2GM / . This will determine 



the constant in the expression for p, Eqn. ()53|) . so that 
for example for a polytrope with F = 5/3, one obtains: 



2GM \ 
^ 5Krout / 



rout 



rout 



3/2 



(56) 



where k = {3/5){7T^/3y/^m~'^^'^ « 2.3 GeV^s/^ ^ 
Given this expression for the density, the effect on the 
tidal acceleration due to the singularity can now be com- 
pletely determined for a given f{R) model. Fig. [6] shows 
the size of i^'sing for f{R) = R — p'^/R, compared to the 
size of the radial tidal acceleration in the Schwarzschild 
solution, Kgb. — —2GM/r^. For typical neutron star 
parameters with M « 2Mq and rout ~ 10 km, one sees 
that the tidal acceleration due to the singularity becomes 
equal to the Schwarzschild value only at a distance ~ 0.3 
fermi from the surface (which still corresponds to a tiny 
acceleration w 3 x 10~® m/s^ over the distance of one 
fermi). It should be clear that extending the validity of 
the polytropic EOS and the mathematical definition of a 
singular surface to these dimensions is completely unreal- 
istic. We conclude that no generic qualitative restrictions 
can be put on Palatini f{R) gravity based on the diver- 
gence of p" for polytropic EOS's at the boundary. 
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Ar / fermi 



FIG. 6: Shown is the logarithm of the ratio Ksing/ Kqk 
for a polytropic star with F = 5/3 in a Palatini f[R) — 
R — j-L'''/R model, as a function of Ar = rout — r, where 
Kgr = —2GM/r^ is the coefficient of radial tidal acceler- 
ation in the Schwarzschild solution. The dashed curve shows 
the corresponding density in units lO^^pcrit, which goes to 
zero at the sharp edge of the star Ar = 0. We have chosen 
M = 2Mq and rout = 10 km. 



As was already pointed out in Ref. [lu\, one does ex- 
pect that quantitative constraints should arise from com- 
pact stars for the Palatini f{R) models where F is an 
increasing function of R. This expectation seems to be 
at least partially realized by the results in Ref. |23j] on 
an f{R) = R + eR^ model. However, the slight deviation 
from the GR metric found to occur inside the star might 
not provide a striking enough an observable difference to 
rule out even these types of models. We plan to return 
to this issue in a future publication. 

Let us finally comment on another recent claim in 
Ref. [1^, according to which Palatini f{R) gravity would 
always lead to infinite tidal forces at the boundaries of 
any objects, such as planets, compact stars, or even 
atoms. On part, the arguments of Ref. [1^ (as was also 
observed in Ref. [lB|) were based on an incorrect identi- 
fication of the freely falling local inertial frames. As to 
the behaviour of gravity near sharp boundaries, our ar- 
guments above apply. Atoms in particular should rather 
be thought of as smooth quantum density distributions 
Patom ~ |'/'atom(2^)P, and low energy gravity corrections 
to such a quantity, either in GR or in a Palatini f{R) 
gravity, are utterly negligible. 



IV. SCALAR-TENSOR EQUIVALENCE 

There is a well known equivalence at the classical level 
between f{R) gravity and Jordan-Brans-Dicke (JBD) 
scalar-tensor theory (for a review see [27]]). The general 



JBD action is given by: 

+ S^[g^,,ib]. (57) 

The metric f{R) theory is then obtained with lo = and 
the Palatini theory with uj = —3/2. In both cases the 
field and its potential [/(</)) are related to the f{R) 
parameters via 

(/) EE F{R) ^ R^n{(j3), 

u{^) = ^n{4>)- f{n{4>)). (58) 

The action ((57|) is written in the Jordan frame, defined 
as the frame where the matter fields couple only to the 
metric defining the volume element of the theory. That 
is, the matter action Sra is independent of (j). Note that 
although the field (p has a "wrong" sign for the kinetic 
term in the Palatini case, it does not mean that the field is 
unstable. In the Jordan frame, the gravity sector consists 
of mixing spin-0 and spin-2 fields, and to study stability 
issues one should move to the Einstein frame. 

The Jordan frame is however most convenient for dis- 
cussing physical observations, because the freely falling 
local inertial frames correspond to the coordinate sys- 
tems with locally flat Jordan frame metric. This is so 
because matter follows the geodesies set by the Jordan 
frame metric q„y [T7| . One can in particular show that 
(see e.g. Ref. [2g]) in a massless Jordan-Brans-Dicke the- 
ory the 7PPN parameter is simply related to lo: 

7PPN = '^^^-^-^ ■ (59) 

It was through this result that Chiba [Si] , using the STG 
equivalence, first found that metric f{R) gravity give 
7ppn = 1/2. Turning the argument around, one sees 
that the PPN bound |1 — 7ppn| ^ 10"'' puts a strong 
limit on the Brans-Dicke parameter: lo > 10^. 

Let us now address the PPN bound for the Palatini 
f{R) gravity in the STG context. Eqn. ([59]) cannot be 
used as such and the limit co — > —3/2 turns out to be 
very delicate. It is instructive to rewrite the theory in 
the Einstein frame by introducing a new metric variable: 

hf^„ = (pg^^ . (60) 

When written in terms of this metric, the action ([57|) 
becomes: 

+ ^m[0~'V>^]- (61) 

where ip = y372Mpi log </>, V{ip) = M^^U {(j)) / 2(t>^ , and 
Mpi = (87rG)"^/^ Note that the actions ^ and ((6T1) 
are completely equivalent; the conformal scaling of the 
metric (|60p is nothing but a convenient change of vari- 
ables to a frame where the spin-0 and spin-2 degrees of 
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FIG. 7: Shown is the effective potential Vefr(<y3) (in units 
H^Mpi) as a function of ip/Mpi = ,/372 log(F/Fmin) for the 
model f{R) ^ R~ ^j.'^/R + {a/2fi'^)R'^. The solid line cor- 
responds to a = 0.8 with p = 5pcrit and the dashed line to 
a = 0.01 with p — O.lpcrit- 



freedom of the gravity sector decouple. There has been 
a surprising amount of confusion related to this issue, 
including debates as to which of the two frames should 
be considered the "physical" one [1^. Of course, neither 
frame (or any other of an infinity of possible choices for 
frames) is any more physical than the other; some phe- 
nomena are easier to interpret in the Jordan frame while 
others are best worked out using the Einstein frame. For 
a clear discussion of these issues see e.g. Refs. [28l . [30| . 
For example, the stability of the JBD theory is obscured 
by the spin mixing in the Jordan frame, but becomes 
obvious in the Einstein frame, where the kinetic terms 
in the action (|61[) are manifestly positive definite given 
uj > —3/2. For Palatini f{R) gravity, the kinetic term 
vanishes in the Einstein frame, and so these models cor- 
respond to a boundary between stable and unstable JBD 
models. 

The equation of motion for ip m a, generic JBD theory 
follows from the action (ICTI): 



extremum of VcS- As is seen in Fig. this extremum 
can correspond to either a minimum or a maximum of the 
potential, depending on the model parameters. For the 
model f{R) — R — IJL^ / R in particular, this extremum is a 
maximum. In the metric version of the ,f{R) theory {uj — 
0), the potential remains the same, but the field is no 
longer constrained to stay at the extremum value. This is 
obviously the origin of the Dolgov-Kawasaki instability in 
the STG language and it affects all models with V^^ < 0. 

Now, consider a JBD theory in the limit uj —3/2. In 
this limit the solution for the trace equation ([55)1 becomes 
an infinitely strong attractor in the class of static solu- 
tions for any model with V^g < 0. However, at the same 
time this solution becomes arbitrarily unstable in time. 
Indeed, perturbing the field (f> around the static solution, 
(f){r) — > 0(r) -I- d(t){t, r), it is straightforward to show that 
a Fourier mode 5(j)k of the linearized perturbation obeys 
the equation 



2uj + 3 



' en 



(64) 



In the limit lo — s- —3/2, all modes become unstable at the 
same time as their characteristic decay time goes to zero 
for V^g < 0. This argument shows that it is not possible 
to obtain a Palatini theory with an effective potential 
V^g < 0, as a limiting case of a JBD theory with lu — > 
—3/2. This is true in particular for the f{R) — R — fi'^/R 
model. For models with V^g > 0, the continuous limit 
does exist however. Nevertheless, a theory with sufficient 
stability is not necessarily acceptable as a Palatini f{R) 
gravity model. For example, stabilizing the vacuum of 
a JBD model with a potential derived from the model 
f{R) = R- fj,^/R + (a/2/i2)i?2^ sets the bound a > 1. 
However, we have seen that in the Palatini limit this 
leads to an enormous value of the function F (see section 
III C ip . giving rise to an unacceptable exterior metric. 
Even though Palatini f{R) theories are obtained, and 
can be considered independently of JBD theories, the fact 
that a smooth limit does not exist for V^g < models is 
troubling and may indicate some fundamental problem 
with such theories. 



2uj 



p 



g-Vs/Sv'/A^pi = , (62) 



PI 

where we have neglected pressure, and the exponential 
factor appears since p is the Jordan frame energy den- 
sity. Now, for Palatini f{R) gravity the coefficient in 
front of the kinetic term is strictly zero, and the result- 
ing equation 



^/6A^J 



/8/3(/>/J\fpi 



0, 



(63) 



PI 



is easily seen to be equivalent with the trace equation 
([37]) . The Palatini field is then completely fixed by the 
potential V{if!) (i.e. the form of the function f{R)) and 
the matter density: the field is constrained to sit at the 



V. CONCLUSIONS 

In this paper we have studied the interior and exterior 
spacetimes of stars in f{R) gravity theories. We started 
by deriving the generalized Tolman-Oppenheimer-Volkoff 
equations for spherical hydrostatic equilibrium, for both 
metric and Palatini versions of the theory. These equa- 
tions were solved analytically in the Newtonian limit, 
and the results were supplemented with numerical cal- 
culations. Our analysis showed that metric f{R) theo- 
ries, for the major part of the parameter space, predict a 
PPN parameter 7ppN = 1/2 in contrast with Solar Sys- 
tem observations. These solutions were independent of 
finite matter configurations (dark matter) surrounding 
the star. 
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We also found a class of solutions corresponding to 
boundary values where Fq — 1 < 10~^, for which 7ppn ~ 
1 in metric f{R) models. However, these solution are 
unstable against decay in time and are hence unphysical. 
Our analytic results were derived mainly in the context of 
the specific f{R) = i? — fi^/R model, but our numerical 
analysis supports the conclusion that the metric f{R) 
models (beyond the Einstein-Hilbert limit f{R) = R — 
2A) either predict 7ppn = 1/2, or are unstable. However, 
the Palatini versions of these theories were found to be 
easily compatible with the Solar System constraints. In 
addition, we showed that compact stars with polytropic 
EOS's are consistent with Palatini f{R) gravity. Thus 
there appears to be no support for a recent claim for a 
"no-go theorem" for Palatini f{R) gravity, based on the 
use of a polytropic EOS [15J . 

We finally considered the equivalence between f{R) 
gravity and scalar-tensor theory. While both the metric 
and the Palatini versions of the theory can formally be 



viewed as a particular case of Jordan-Brans-Dicke theory, 
understanding Palatini /(i?) gravity as a limiting case of 
a Jordan-Brans-Dicke theory with uj —3/2, is not pos- 
sible at least for models with V^^ < 0. This is somewhat 
troubling and may indicate some fundamental problem 
with such theories. 
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